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ABSTRACT 

We use a semi-analytical model for the substructure of dark matter haloes to assess the 
too-big-to-fail (TBTF) problem. The model accurately reproduces the average subhalo 
mass and velocity functions, as well as their halo-to-halo variance, in TV-body simu¬ 
lations. We construct thousands of realizations of Milky Way (MW) size host haloes, 
allowing us to investigate the TBTF problem with unprecedented statistical power. We 
examine the dependence on host halo mass and cosmology, and explicitly demonstrate 
that a reliable assessment of TBTF requires large samples of hundreds of host haloes. 
We argue that previous statistics used to address TBTF suffer from the look-elsewhere 
effect and/or disregard certain aspects of the data on the MW satellite population. 
We devise a new statistic that is not hampered by these shortcomings, and, using only 
data on the 9 known MW satellite galaxies with V max > 15kms _1 , demonstrate that 
1.4jli 1 % of MW-size host haloes have a subhalo population in statistical agreement 
with that of the MW. However, when using data on the MW satellite galaxies down 
to Vmax = 8 kms -1 , this MW consistent fraction plummets to < 5 x 10 -4 (at 68% 
CL). Hence, if it turns out that the inventory of MW satellite galaxies is complete 
down to 8 kms^ 1 , then the maximum circular velocities of MW satellites are utterly 
inconsistent with ACDM predictions, unless baryonic effects can drastically increase 
the spread in V max values of satellite galaxies compared to that of their subhaloes. 

Key words: methods: analytical — methods: statistical — galaxies: haloes — cos¬ 
mology: dark matter — Galaxy: halo 


1 INTRODUCTION 

The A+cold dark matter (ACDM) paradigm of structure 
formation has been remarkably successful in explaining 
cosmic structure covering a wide range in redshift and 
scale. However, on small scales, at low redshifts, a num¬ 
ber of potential problems have been identified regarding 
the abundance and/or structural properties of dark mat¬ 
ter (sub)haloes and their associated (satellite) galaxies. In 
the order in which they have been introduced, these are the 
‘cusp-core’ problem, according to which the cuspy central 
density profiles of CDM haloes predicted by dark-matter- 
only simulations are inconsistent with observed rotation 
curves (e.g., Flores & Primack 1994; Moore 1994; Kuzio de 
Naray, McGaugh & de Blok 2008; Trachternach et al. 2008; 
de Blok 2010; Oh et al. 2011; but see also van den Bosch & 
Swaters 2001 and Dutton et al. 2005), the ‘missing-satellite’ 
problem, highlighting the large discrepancy between the pre¬ 
dicted number of CDM sublialoes per host halo and the 
much smaller number of satellite galaxies detected around 
host galaxies such as the Milky Way and M31 (e.g., Klypin 
et al. 1999; Moore et al. 1999), and the ‘too big to fail’ 
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(TBTF) problem, which basically refers to the overabun¬ 
dance of massive, dense subhaloes predicted by CDM com¬ 
pared to the observed number of relatively luminous satellite 
galaxies of the Milky Way or the Local Group (e.g., Boylan- 
Kolchin, Bullock & Kaplinghat 2011, 2012; Martinez 2013; 
Garrison-Kimmel et al. 2014a, Tollerud, Boylan-Kolchin & 
Bullock 2014). 

In this paper we focus on the TBTF problem (hereafter 
simply ‘TBTF’), which is generally considered the most diffi¬ 
cult to reconcile with ACDM, and which has spurred a frenzy 
of papers offering possible solutions and/or advocating mod¬ 
ifications of the standard paradigm. This includes sugges¬ 
tions to change the nature of the dark matter from ‘cold’ to 
either ‘warm’ or ‘self-interacting’ (e.g., Maccfo & Fontanot 
2010; Vogelsberger, Zavala & Loeb 2012; Lovell et al. 2012; 
Anderhalden et al. 2013; Rocha et al. 2013; Shao et al. 2013; 
Polisensky & Ricotti 2014), relatively small changes in the 
normalization, ag, and/or spectral index, n a , of the initial 
power spectrum (e.g., Polisensky & Ricotti 2014), a highly 
stochastic star formation efficiency for galactic subhaloes, so 
that a fraction of the more massive subhaloes remain dark 
(e.g., Kuhlen, Madau & Krumholz 2013; Rodriguez-Puebla, 
Avila-Reese & Drory 2013a,b), lowering the mass of the MW 
host halo to ~ 10 11,8 /i _1 Mq (Di Cintio et al. 2011; Wang et 
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al. 2012; Vera-Ciro et al. 2013), and enhanced tidal (impul¬ 
sive) heating of satellite galaxies due to the stellar disk of 
the Milky Way (Zolotov et al. 2012; Brooks & Zolotov 2014; 
Arraki et al. 2014). 

Despite all this interest in TBTF, we still lack consensus 
regarding both its formulation and its severity. In fact, the 
TBTF problem has been formulated in a few subtly different 
ways. For example, Wang et al. (2012) express it as a prob¬ 
lem of missing massive satellites; high-resolution A-body 
simulations of Milky-Way (MW) size host haloes reveal of 
the order of 10 subhaloes per host with a maximum circular 
velocity Vmax (J, 25kms _1 (e.g., Strigari et al. 2007; Madau, 
Diemand & Kuhlen 2008; Boylan-Kolchin et al. 2011). Yet, 
only two MW satellite galaxies (the large and small Magel¬ 
lanic clouds, hereafter LMC and SMC) are believed to have 
associated subhaloes that meet this criterion. We shall refer 
to this as the ‘massive subhaloes formulation’ of TBTF. A 
related expression of TBTF is that the MW reveals a gap in 
the Knax-distribution of its satellite galaxies. There are no 
satellite galaxies known with 25kms _1 < V max ^ 55 km s -1 , 
something that is not reproduced in numerical simulations 
(e.g., Boylan-Kolchin et al. 2012, Cautun et al. 2014b). We 
shall refer to this as the ‘gap formulation’ of TBTF. Finally, 
in what we call the ‘density formulation’ of TBTF, it is ar¬ 
gued that the (central) densities of dark matter subhaloes 
are too high compared to the central densities in satel¬ 
lite galaxies as inferred from their kinematics (e.g., Boylan- 
Kolchin et al. 2011; Purcell & Zentner 2012). In this paper 
we will compare all three formulations, and point out how 
and where they differ. In particular, we will point out how 
some of the solutions listed above may solve TBTF in one 
formulation, but not in the other(s). 

As to the severity of TBTF, the lack of consensus largely 
owes to poor statistics. On the observational side, since 
TBTF relates to low-luminosity satellite galaxies, which can 
only be observed in the local neighborhood, the observa¬ 
tional ‘evidence’ for a TBTF problem is limited to the MW 
and M31. Ideally, one would like to test how common this 
problem is in other host haloes, spanning a range in halo 
masses and environments. On the theoretical side, the TBTF 
problem was originally identified using a sample of only 6 
simulated MW-size host haloes from the Aquarius project 
(Springel et al. 2008). More recently, the sample size of high 
resolution MW-size haloes has grown to the order of 100, 
with the ELVIS (Garrison-Kimmel et al. 2014a) and cl25- 
2048 (Mao, Williamson & Wechsler 2015) suites each adding 
of order 50. However, as we demonstrate in this paper, accu¬ 
rately capturing the halo-to-halo variance, and characteriz¬ 
ing the severity of TBTF at percent level accuracy, requires 
of order a thousand realizations of MW-size host haloes 
with subhaloes resolved down to Vmax ~ 10 kins' 1 . This 
is immensely challenging, even for the largest supercomput¬ 
ers available to date (see discussion in van den Bosch et 
al. 2014). 

Because of these statistical drawbacks, some authors 
have resorted to semi-analytical models (Purcell & Zenter 
2012) or empirical extrapolations of large simulations (Cau¬ 
tun et al. 2014a,b) to generate large samples of well-resolved 
MW-size haloes. In particular, Purcell & Zenter (2012; PZ12 
hereafter) use a model developed by Zentner et al. (2005) 
to generate thousands of MW-size haloes, and claim that at 
least ~ 10% of MW-size haloes are free from the TBTF prob¬ 


lem (using the density formulation). However, this study 
suffers from three drawbacks. First of all, the halo merger 
trees they use, and which are the backbone of their analyti¬ 
cal model, are constructed using the method of Somerville & 
Kolatt (1999), which has been shown to be significantly and 
systematically biased (Zhang, Fakhouri & Ma 2008; Jiang 
& van den Bosch 2014a). Second, their recipe for computing 
Vjnax of subhaloes is oversimplified, resulting in unrealistic 
density estimates. And finally, they only address TBTF in 
the density formulation, and it therefore remains to be seen 
what their model predicts regarding the other two formula¬ 
tions. 

In this study, we revisit the TBTF problem using a 
semi-analytical model that we developed and tested in Jiang 
& van den Bosch (2014b; hereafter Paper I) and van den 
Bosch & Jiang (2014; hereafter Paper II). The model uses 
accurate halo merger trees, and well-calibrated (yet simple) 
recipes for the tidal stripping and disruption of subhaloes. 
We demonstrate that the model accurately reproduces the 
halo-to-halo variance of subhalo properties found in numer¬ 
ical simulations, and use it to gauge the severity of TBTF 
in all three formulations discussed above. We examine the 
dependence on host halo mass and cosmology, and explicitly 
demonstrate that numerical simulation suites with of order 
50 host haloes can only assess TBTF statistics at ~ 10 per¬ 
cent accuracy. Finally, we show that all three formulations 
of TBTF either are hampered by the ” look-elsewhere effect” 
or disregard certain aspects of the data. In order to remedy 
these shortcomings, we devise a new statistic that treats all 
data on equal footing, without the need to pre-identify spe¬ 
cific Vmax scales in the data. Application of this new statis¬ 
tic to the existing data on the MW satellite galaxies makes 
it clear that TBTF is predominantly a statement about the 
Vmax distribution of MW satellites being much broader than 
predicted by the ACDM paradigm. 

The paper is organized as follows. In §2 we briefly de¬ 
scribe the model and demonstrate that it generates subhalo 
populations that are indistinguishable from those of large 
N-body simulations. §3 employs thousands of model real¬ 
izations to evaluate what fraction of MW-size haloes has 
subhalo populations consistent with the satellite properties 
of the Milky Way, using all three formulations of the TBTF 
problem. In §4 we discuss the severity of TBTF in light of 
our new statistic, and we summarize our findings in §5. 


2 SEMI-ANALYTICAL MODEL 

The basis of this work is a semi-analytical model that we 
developed in Paper I. The model is designed to generate 
subhalo populations for a target host halo of any mass, in 
any reasonable ACDM cosmology. This section describes the 
model, and shows that it yields subhalo statistics in excel¬ 
lent agreement with state-of-the-art numerical simulations. 
However, we start with a brief introduction of halo basics, 
outlining a number of definitions and notations. 

2.1 Halo Basics 

Throughout this paper, dark matter haloes at redshift « are 
defined as spherical systems with a virial radius r v ir inside 
of which the average density is equal to A v i r (z)p cr it(z). Here 
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pait{z) = 3H 2 (z)/8tvG is the critical density, and A v i r (z) is 
given by 

A vir (z) = 18tt 2 + 82* - 39* 2 (1) 


with * = fi m (z) — 1 (Bryan & Norman 1998). Haloes whose 
center falls within the virial radius of another halo are called 
subhaloes, while haloes that are not subhaloes are called 
host haloes. The (virial) mass of a host halo is defined as 
the mass within the virial radius r vlI and indicated by M. 
The mass of a subhalo is defined as the fraction of its mass 
at accretion (i.e., when it transits from being a host halo to 
a subhalo) that remains bound, and is indicated by m (see 
Paper II for an in-depth discussion of subtleties associated 
with mass definitions). 

Throughout we assume that host haloes have a NFW 
density profile (Navarro, Frenk & White 1997), characterized 
by a concentration parameter c = r v ii/r a , with r s the NFW 
scale radius, and a maximum circular velocity of 


Fmax = 0.465Kir 

Here 


ln(l + c) — c/( 1 + c) ’ 


( 2 ) 


Kir = = 159 - 43 kms -1 ^ f 0l2 ^_im q 

1 1/6 


M 


1/3 


H(z) 1 

1/3 

^vir (-^) 

Ho 


178 


(3) 


is the virial velocity, and the radius, r ma x, at which the cir¬ 
cular velocity reaches its maximum value, is given by 


^max — 2.163 r s . 


(4) 


2.2 Model Description 


The semi-analytical model that we use to construct real¬ 
izations of dark matter subhalo populations is a strongly 
improved and modified version of the model introduced in 
van den Bosch, Tormen & Giocoli (2005; hereafter vdB05). 
It treats subhalo mass stripping in an orbit-averaged sense, 
allowing the construction of thousands of realizations in a 
manner of minutes. This section gives a brief overview of 
the model ingredients. A more detailed description can be 
found in Paper I. 

The backbone of our model, as for any other semi- 
analytical model for the substructure of dark matter haloes 
(e.g., Taylor & Babul 2001, 2004, 2005a,b; Benson et 
al. 2002; Taffoni et al. 2003; Oguri & Lee 2004; Zentner 
& Bullock 2003; Pcnarrubia & Benson 2005; Zentner et 
al. 2005; vdB05; Gan et al. 2010; Yang et al. 2011; Pur¬ 
cell & Zentner 2012), is halo merger trees. These describe 
the hierarchical mass assembly of dark matter haloes, and 
therefore yield the masses, m ac c, and redshifts, z acc , at which 
the dark matter subhaloes are accreted into their hosts. We 
also use the merger trees to compute the concentration pa¬ 
rameter, c a cc, of the subhalo at accretion. For this we use 
the fact that halo concentration is tightly correlated with 
its assembly history (e.g., Wechsler et al. 2002; Ludlow et 
al. 2013). In particular, we compute c aC e using the model of 
Zhao et al. (2009), according to which 


c(t) = 4.0 



t 


L3.75to.o4 J 


1/8 


(5) 


Here to .04 is the time at which the main progenitor of the 
halo in question assembled 4 percent of its mass at time 
t, which we compute from our merger tree by tracing the 
assembly history of the subhalo back in time. Using eqs. (2)- 
(4), combined with m aC c, « ac c and c acc , we also compute U max 
and r max of the subhalo at accretion, which we indicate by 
Kcc and r aC c, respectively. 

After accretion, a subhalo orbits its host halo, subject 
to tidal stripping, impulsive heating and dynamical friction. 
Following vdB05 we model the tidal stripping in an orbit- 
averaged sense, where the average is taken over all orbital 
energies, eccentricities and phases. Using a simple toy model, 
in which it is assumed that over a radial orbital period the 
subhalo is stripped of all mass outside of the subhalo’s tidal 
radius at the orbit’s pericentric distance from the center of 
the host halo (see Paper I for details), one predicts an orbit- 
averaged mass loss rate that is well described by 



with m, M, and Td yn the instantaneous subhalo mass, host 
halo mass, and host halo’s dynamical time, respectively. The 
same functional form was also used by Giocoli, Tormen & 
van den Bosch (2008) to fit the orbit averaged mass loss 
rates of dark matter subhaloes in numerical IV-body simu¬ 
lations, which yielded A = 1.54^q' 3 2 and ( = 0.07 ± 0.03. 
Our toy model predicts that ( ~ 0.04 and that A follows a 
log-normal distribution with median A — 0.81 and scatter 
a\ og A — 0.17. This scatter in A is due to the variance in or¬ 
bital properties (energy and angular momentum) and halo 
concentrations (of both the host and the subhalo). Given 
the uncertainties in the parameters derived from the sim¬ 
ulations, and the oversimplifications of the toy model, we 
treat A and (j as free parameters. As detailed in Paper I, we 
tune A and ( so as to reproduce the subhalo mass functions 
in the Bolshoi simulation (see §2.3 below). This results in 
A = 0.86 and £ = 0.07, close to the values suggested by 
the toy model, and in good agreement with the simulation 
results of Giocoli et al. (2008). The scatter we keep fixed 
at <7i 0 g a = 0.17. For each individual subhalo in our model, 
we randomly draw a value for A from the log-normal, and 
evolve its mass using Eq. (6). 

During its pericentric passage, a subhalo experiences 
impulsive heating which increases the kinetic energy of its 
constituent particles. Depending on the amount of energy 
transferred, the subhalo either expands (‘puffs up’) while re¬ 
establishing virial equilibrium, or is tidally disrupted (typi¬ 
cally if the kinetic energy injected exceeds the gravitational 
binding energy of the subhalo). In our model, a subhalo is 
disrupted (i.e., removed from the inventory) if its mass m(t) 
drops below a critical mass 

m-dis — niacc(< /dis*s, a cc), (7) 

with m ac c(< r) the mass enclosed within radius r at ac¬ 
cretion, and r s , a cc the NFW scale radius of the subhalo at 
accretion. The dependence on r s , a cc assures that subhaloes 
that are more concentrated at accretion are more resistant 
to disruption. Using idealized IV-body simulations, Hayashi 
et al. (2003) found that /di s — 2, while both Taylor & Babul 
(2004) and Zentner et al. (2005) used Eq. (7) to model 
tidal disruption in their semi-analytical models, but with 
wildly differing values of /dis = 0.1 and /dis = 1-0, respec- 
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tively. As detailed in Paper I, tidal disruption in numerical 
./V-body simulation occurs for values of fdis that follow a 
broad (roughly log-normal) distribution with a median of 
fdis ~ 1.5 and standard deviation criog(/ di ) ~ 0.55. Based 
on these findings, for each subhalo we draw a value of fdis 
from a similar log-normal distribution, and disrupt the sub¬ 
halo whenever its mass drops below mdi S . By calibrating 
the model to reproduce the distribution of retained mass 
fractions, m/m ac c, of the surviving subhaloes in the Bol¬ 
shoi simulation, we find that we need to introduce a weak 
mass dependence in fdis, and we end up modeling the fdis 
distribution as a log-normal with median 


fdis = 1.5 


1 + 0.8 




(8) 


and standard deviation cri og (y ) = 0.55. 

The model for mass stripping and tidal disruption de¬ 
scribed above regulates the mass evolution of the subhaloes. 
However, in order to address TBTF, we also need to model 
the structure (i.e., density distribution) of the individual 
subhaloes. In particular, we need to be able to predict V max 
and r max . Based on the idealized simulations of Peiiarrubia 
et al. (2008, 2010), we model y max and r max of the subhaloes 
using 


Pm ax 
Vacc 


a: 0 ' 3 

1.32-;-rrr—r 

(1 + a:) 0 - 4 


T max 
T acc 


I81 (l + .r)“°- 3 ’ (9) 


with x = m/rria.ee', i.e., the evolution of V ma x and r max of 
subhaloes depends only on the amount of mass loss, but not 
on how that mass was lost (see also Hayashi et al. 2003). As 
shown in Paper I, Eq. (9) is in good agreement with results 
from various numerical simulations. 


2.3 N-body Simulations 

In this paper we use several numerical N-body sim¬ 
ulations for comparison with our model predictions. 
The first is the Bolshoi simulation (Klypin, Trujillo- 
Gomez & Primack 2011), which follows the evolu¬ 
tion of 2048 3 dark matter particles in a box of size 
250 ft -1 Mpc using the Adaptive Refinement Tree (ART) 
code (Kravtsov, Klypin & Khokhlov 1997) in a flat ACDM 
cosmology with parameters (fi m ,o, Da,o, fib.o, h, as, n B ) = 
(0.27, 0.73, 0.047, 0.7, 0.82, 0.95) (hereafter ‘Bolshoi cosmol¬ 
ogy’)- With a particle mass m p = 1.4 x 10 8 /!+ 1 Mq, Bol¬ 
shoi resolves haloes down to ~ 10 10 /i -1 Mq (correspond¬ 
ing to Vmax ~ 50 knis -1 ). We use the publicly available 
halo catalogs^ obtained using the phase-space halo finder 
R0CKSTAR (Behroozi et al. 2013a,b). R0CKSTAR haloes are de¬ 
fined as spheres with an average density of A v i r p C rit, in line 
with the halo definition used throughout this paper. 

In addition to the Bolshoi simulation, we also com¬ 
pare our results to the ELVIS suite of zoom-in sim¬ 
ulations of 48 MW-size dark matter haloes (Garrison- 
Kimmel et al. 2014a). These cover the mass range 
11.85 < log[M/( A -1 Mq)] < 12.31, comparable to the 
range of masses quoted for the Milky Way in the lit¬ 
erature. The haloes have been selected from medium- 
resolution cosmological volumes of box size 50 h -1 Mpc 


and re-simulated with progressively higher resolution up 
to m p = 1.35 x 10 5 h~ 1 M©. The ELVIS suite adopts a 
flat ACDM cosmology with (fl mj o, Da,o, I2b,o> h, as, n s ) = 
(0.266,0.734,0.045,0.7,0.801,0.963) (hereafter ‘WMAP7 
cosmology), which are the parameters that best fit the 7- 
year data release of the Wilkinson Microwave Anisotropy 
Probe (Larson et al. 2011). The halo catalogs of the ELVIS 
suite are publicly available'!', and contain subhaloes (identi¬ 
fied with R0CKSTAR , and using the same halo definition as for 
Bolshoi), down to the resolution limit of Vna,x = 8 kms -1 . 


2.4 Comparison with Simulations 

In this section, we first summarize the results of Paper I & II 
by showing that the model accurately matches the average 
subhalo mass and Vn ax functions extracted from A+body 
simulations, and then demonstrate that the same model, 
without any modifications, also accurately reproduces the 
halo-to-halo variance of subhalo statistics. 

The upper panels of Fig. 1 compare the model 
predictions for the average subhalo mass function, 
dAf/dlog(m/JV/o), and the average subhalo velocity func¬ 
tion, dAl/dlog(Lmcix/Vvir,o) (solid and long-dashed curves), 
with the results from the Bolshoi simulation (filled and open 
circles). The latter are averaged over a total of 281 host 
haloes with mass Mo = io 14 ' 25±0 ' 25 h~ x Mg, while the model 
predictions are obtained averaging over 2000 realizations of 
host haloes with mass Mo = 10 14 ' 25 h -1 Mg, adopting the 
Bolshoi cosmology. Solid lines and filled circles indicate the 
abundances of the surviving subhaloes as a function of their 
present-day mass and Vmax, while dashed lines and open 
circles show the abundances of the same subhaloes, but as 
function of their mass and Vmax at accretion. In both cases, 
the model predictions are in excellent agreement with the 
Bolshoi results. Note that the abundance functions at ac¬ 
cretion are not to be confused with the unevolved subhalo 
mass and/or V max functions, and which are indicated by 
the short-dashed curves. The latter include all subhaloes 
that have ever been accreted onto the main progenitor of 
the host halo, and includes those sublialoes that have been 
disrupted since. Note that the unevolved subhalo mass and 
velocity functions are substantially higher than those for the 
surviving population, indicating that subhalo disruption is 
extremely efficient and important. 

The lower, left-hand panel of Fig. 1 plots the concentra¬ 
tions at accretion of the present-day, surviving subhaloes as 
a function of their mass at accretion. The model prediction 
is again in excellent agreement with Bolshoi, indicating that 
our model accurately reproduces the structural properties 
of subhaloes at infall. Finally, the lower, right-hand panel of 
Fig. 1 compares the distributions of the retained mass frac¬ 
tion, m/m a cc, of the present-day, surviving subhaloes. Note 
that very few subhaloes lose more than 90% of their initial 
mass, which is another manifestation of efficient subhalo dis¬ 
ruption. Here the good agreement between the model and 
Bolshoi lays the foundation for the accurate model predic¬ 
tion of the parameters that are important for TBTF, V max 
and r max , whose evolution is purely a function of m/m ac c 
(see Eq. [9]) 


t http://www.slac.stanford.edu/~behroozi/Bolshoi_Catalogs/ 


+ http://localgroup.ps.uci.edu/elvis/data.html 
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Figure 1. The upper panels compare the average subhalo mass (left) and velocity (right) functions of the Bolshoi simulation (circles) and 
our model (lines). The line- and symbol- styles differentiate the results for the present-day surviving subhaloes (‘present’), the surviving 
subhaloes at accretion (‘accretion’), and all subhaloes ever accreted (‘unevolved’). The lower, left-hand panel plots the concentrations 
of the surviving subhaloes at accretion, Cacc, as a function of their masses at accretion (normalized to the present-day host halo mass). 
Symbols (Bolshoi) and solid curve (model) represent the median relations, while error bars and hatched region indicate the 68% confidence 
intervals. Finally, the lower, right-hand panel shows the distributions of the retained-mass fractions of the surviving subhaloes, with 
downward triangles indicating the medians. 


The set of comparisons shown in Fig. 1 are the key 
diagnostics that we used in Paper I to calibrate our model. 
Note, though, that the comparison is for host haloes with 
Mo ~ 10 14 ' 25 r'Mg, far from the mass scale of interest 
for the TBTF problem. The reason for showing this mass 
scale is that, in the Bolshoi simulation, it probes a large 
dynamic range in subhalo mass (down to m = 10~ 4 Mo) and 
a sufficiently large sample size within the simulation box to 
be able to compute a reliable average mass function. As we 
have demonstrated in detail in Papers I and II, the same 
model is equally successful at other mass scales, as far as 
they have been probed by simulations. 

Fig. 2 plots the individual , cumulative subhalo velocity 


functions, N(> Vmax), for 1986 host haloes with mass Mo = 
1q12.io±o.oi /j-!y 0 (left-hand panel) and 441 host haloes 
with mass Mo = io 13 - 50±0 05 (right-hand panel) in 

the Bolshoi simulation, and compares them with 2000 model 
realizations of halo mass Mo = 10 12 ' 10 /i _1 Mq, and 500 
model realizations of halo mass Mo = 10 13 ’ 50 re¬ 

spectively. The Bolshoi results are plotted down to Vmax = 
50 kms , which corresponds to roughly 250 particles, the 
minimum number of particles required to resolve haloes well 
enough for a reliable estimate of V m ax (see Paper II). For the 
model realizations, we trace subhaloes down to a mass of 
10 _5 Mo.Inthe case of host haloes with Mo = 1O 121O /i _1 M0 
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v max [kms -1 ] V max [kms *] 

Figure 2. Cumulative subhalo velocity functions for individual host haloes. Upper, left-hand panel: 2000 model realizations of haloes 
with Mo = 10 1210 fe-'Mg (grey) and 1986 Bolshoi haloes with Mo = 10 12 - 10±0 01 fe -1 Mg (red). Upper, right-hand panel: 500 model 
realizations of haloes with Mo = 10 13 - 50 fe -1 Mg (grey) and 441 Bolshoi haloes with Mo = lo 13 - 50±005 /i _ 1 Mq (red). The vertical, 
dotted lines mark the Bolshoi resolution limit at Knax = 50 kms -1 . The thick, grey line indicates the model’s median Vmax at given 
N[> y max ), with the dash-dotted lines bracketing the 68% and 95% intervals, as indicated. The thick, red curve (upper, right-hand) is 
the median relation for the Bolshoi simulation. The thick, green curve (upper, left-hand) is the best-fit median relation for the ELVIS 
simulation, as given by Garrison-Kimmel et al. (2014a). Bottom panels: Vmax distribution of the N* ’ subhalo rank-ordered in Vmax- 


this roughly translates to Vmax = 10 kms 1 , which is more 
than sufficient for a detailed assessment of TBTF. 

The model predictions manifest excellent agreement 
with the Bolshoi results for both mass scales. In the case 
of Mo ~ 10 13,50 /j^Mq, for which the median in Bolshoi 
can be measured over an appreciable range in Vinax, the en¬ 
semble average of the model predictions (indicated by the 
thick gray line) agrees almost perfectly with that of the Bol¬ 
shoi simulation data (solid red line). In terms of the halo- 
to-halo variance, the model and Bolshoi results are almost 
indistinguishable down to the resolution limit of the sim¬ 
ulation. This is nicely illustrated in the bottom panels of 
Fig. 2, which show the Vmax distributions at several values 
of N(> Knax), as indicated. 

The average, cumulative subhalo velocity functions are 
usually described by a simple power law, N(> V’max) oc Emax 
(at least for N ■> 2). The Bolshoi simulation, against which 
our model is calibrated, yields a slope of a ~ 2.9 with, as far 


as we can tell given the limited mass resolution, no depen¬ 
dence on host halo mass (see Papers I & II). This is in excel¬ 
lent agreement with our model predictions, which also yields 
a « 2.9 without any significant dependence on host halo 
mass. The ELVIS suite, however, seems to predict a slope 
that is significantly steeper, with a = 3.3 (Garrison-Kimmel 
et al. 2014a). It is unclear what the cause is of this discrep¬ 
ancy. We emphasize that most other simulations all suggest 
that a <, 3 (e.g., Klypin et al. 1999; Reed et al. 2005; Madau 
et al. 2008; Diemand et al. 2008; Wu et al. 2013). The only 
simulations for which similarly large values for a have been 
reported are the Millennium-II simulation (Boylan-Kolcliin 
et al. 2009) and the suite of Aquarius simulations (Springel 
et al. 2008). Both of these were analyzed with the halo 
finder SUBFIND (Springel, White & Hernquist 2001), which 
has been shown to give much steeper subhalo mass and ve¬ 
locity functions compared to other halo finders, including 
R0CKSTAR (see Paper II). Based on these considerations, we 
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Figure 3. Left-hand panel: the fraction of host haloes < IVms massive subhaloes for the 4800 model realizations of ELVIS-size haloes 
(thick, grey line), the actual ELVIS simulation suite (green), and the 100 mock ELVIS suites (thin, grey lines). Realizations with 7 Vms < 2 
are considered as MW-consistent. Middle panel: same as the left -hand panel, but for 10,000 realizations of different halo masses and 
cosmologies, as indicated. Right-hand panel: the probabilities, i.e., fraction out of 10,000 realizations, of having no more than two massive 
subhaloes (grey) and no less than two Magellanic-Cloud analogs (orange) as a function of halo mass. The upper and lower bounds of 
each band correspond to different threshold V m ax values, as indicated, while solid and dashed curves indicate results for the WMAP7 
and Planck cosmologies, respectively. 


trust the (current calibration of the) model for the study of 
MW-size hosts. If future simulations confirm steeper power- 
law slopes at Galactic mass scales, then it means that some 
of our model parameters, such as £ and A, must have some 
mass dependence §, for which we thus far see no indication in 
the Bolshoi simulation that covers 2-3 orders of magnitude 
in host halo mass. 


3 ASSESSING THE TOO BIG TO FAIL 

PROBLEM 

In this section, we use the model introduced above to gauge 
the severity of TBTF in the massive subhaloes formulation 
(§3.1), the gap formulation (§3.2), and the density formula¬ 
tion (§3.3). We use our model to generate 100 realizations 
of the ELVIS suite as follows. Adopting the same WMAP7 
cosmology as used for the ELVIS simulations, for each halo 
mass in the ELVIS suite of 48, we generate 100 model re¬ 
alizations with a mass resolution of m/Mo = 10 -5 . The 
resulting ensemble of 4800 MW-size host haloes (Mo = 
1q12.08±o. 23 /j-ijyfg) j s then split in 100 mock ELVIS suites, 
each with exactly the same distribution of halo masses. We 
use this set below to compare our model predictions with 
those from the ELVIS suite, and to gauge the suite-to-suite 
variation of the various TBTF statistics. In addition, in or¬ 
der to probe the dependence on halo mass and cosmological 
parameters, we also construct ensembles of 10,000 model 
realizations each for several values of Mo, and for a few dif¬ 
ferent cosmologies. 

When comparing model predictions to data, we use 
the Vnax values for the MW and its satellites listed in Ta¬ 
ble 1. These are compiled from Xue et al. (2008), van der 
Marel & Kallivayalil et al. (2014), Kallivayalil et al. (2013), 

§ For example, we find that with f = 0.19 and A = 1.29 the 
model accurately reproduces the ELVIS results. 


Kuhlen (2010), and Boylan-Kolchin et al. (2012). Most of 
the V max values for the dSphs are taken from Kuhlen (2010) 
and Boylan-Kolchin et al. (2012). Both studies used similar 
methodology which we briefly describe in what follows. 

The only kinematic information available for the MW 
dSphs are the line-of-sight velocities of individual stars, 
which constrain the dynamical mass, M(< r), enclosed 
within some radius, r, that is typically much smaller than 
v max (e.g., Strigari et al. 2008; Wolf et al. 2010). In order 
to infer the corresponding Vnax, Kuhlen (2010) and Boylan- 
Kolchin et al. (2012) assign weights to subhaloes in the Via 
Lactea II and Aquarius simulations, respectively, based on 
how closely they match the measured M(< r). The corre¬ 
sponding dSph is then assigned the weighted-average V ma x 
of these subhaloes. A crucial underlying assumption of this 
method is that the subhaloes in the zoom-in simulations 
used are representative of the ACDM subhalo populations 
of MW sized haloes. The validity of this assumption is chal¬ 
lenged by the dramatic halo-to-halo variance of subhalo pop¬ 
ulations. Fortunately, for the few dSphs covered by both 
Kuhlen (2010) and Boylan-Kolchin et al. (2012), the inferred 
Vnax values are mutually consistent within the errors. When 
available, we use the more recent results of Boylan-Kolchin 
et al. (2012), since they are based on a larger sample of six 
host haloes, as compared to only one in the case of Kuhlen 
(2010). 

For the few dwarfs without published Vnax constraints 
(Sgr dwarf, Bootes II, Segue II, Bootes I, and Leo V), we 
use the relation Vnax = 2.2<tlos advocated by Rashkov et 
al. (2012), and the published stellar line-of-sight velocity dis¬ 
persion ctlos (McConnachie 2012 and references therein). 
The resulting Vnax values are indicated in brackets in Ta¬ 
ble 1. Note that we do not include the Canis Major and 
Bootes III stellar overdensities in this list, as we consider 
them already disrupted. We emphasize, though, that includ¬ 
ing them in the inventory of MW satellites has no impact 
on any of our results and/or conclusions. 
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Table 1. Maximum circular velocity of MW and its satellites. 


Object 

Vmax 
[ kms -1 ] 

Reference 

Milky Way 

170.0 ± 15.0 

[i] 

LMC 

91.7 ± 18.8 

[2] 

SMC 

60.0 ±5.0 

[3] 

Sagittarius 

(25.1 ± 1.5) 


Bootes II 

(23.1 ± 16.3) 


Draco 

20.5±« 

[5] 

Ursa Minor 

20.0^22 

[5] 

Fornax 

17.8 ±0.7 

[5] 

Sculptor 

1 7 Q + 2 - 2 

1 ' • -2.0 

[5] 

Leo I 

16 4“*" 2 - 3 
±o -^-2.0 

[5] 

Ursa Major I 

14+ 8 

[4] 

Ursa Major II 

13tS 

[4] 

Leo II 

12.8± 2 'g 

[5] 

Sextans 

H+o.9 

[5] 

Canes Venatici I 

n-slJl 

[5] 

Carina 


[5] 

Canes Venatici II 

H+2 

[4] 

Hercules 

H+ 3 

1 —1.6 

[4] 

Segue I 

io+ 7 

1U —1.6 

[4] 

Coma Berenices 

9 1+ 2 - 9 

y, -0.9 

[4] 

Willman 1 

o o+2.7 

o-^-o.s 

[4] 

Leo V 

(8-i±!:i) 


Segue II 

(7.5±*-|) 


Bootes I 

(5.3j+?) 


Leo IV 

5+0.8 

[4] 


The references for the values of V m ax listed in Column (3) are: 
[1] Xue et al. (2008); [2] van der Marel & Kallivayalil (2014); 
[3] Kallivayalil et al. (2013); [4] Kuhlen (2010); and [5] Boylan- 
Kolchin et al. (2012). Values in brackets are inferred from the 
empirical relation for MW dSphs, V max = 2.2<tlos (Rashkov et 
al. 2012), with <tlos measurements from McConnachie (2012) 
and references therein. 

3.1 The Abundance of Massive Subhaloes 

The TBTF problem was originally expressed as a ten¬ 
sion between the rotation curves of the ~10 most mas¬ 
sive subhaloes in MW-size host haloes and the kinemat¬ 
ics data of the ^10 brightest MW dwarf spheroidals (here¬ 
after dSphs). The MW dSphs all have stellar kinemat¬ 
ics consistent with Vlnax 25 knrs -1 , while the subhaloes 
have Vmax 25kins -1 . Therefore, several studies have used 
the abundance of subhaloes with Vlnax greater than some 
threshold value (typically in the range of 25-30 kms -1 ) 
as a measure of the TBTF severity (e.g., Boylan-Kolchin 
et al. 2011, 2012; Wang et al. 2012; Garrison-Kinrmel et 
al. 2014b). In particular, Garrison-Kimmel et al. (2014b) de¬ 
fine ‘massive failures’ as subhaloes that started out massive 
(Kcc > 30 kms -1 ) and remain massive (Vmax > 25 knrs -1 ) 
to the present day. The argument is that such subhaloes have 
potential wells in which galaxy formation is expected to ‘suc¬ 
ceed’, to the extent that the absence of a satellite galaxy sig¬ 
nals a ‘massive failure’. We refrain from this nomenclature, 
as it leads to confusion when addressing the LMC and SMC; 
instead, we simply refer to subhaloes with 14cc > 30kms -1 
and Vmax > 25 kms -1 as ‘massive subhaloes’. 

The left-hand panel of Fig. 3 plots the cumulative dis¬ 


tribution of the number of massive subhaloes, Ams, per 
host halo, in each of the 100 mock Elvis suites (thin, gray 
lines). For comparison, the thick, green line shows the results 
from the actual ELVIS suite. Typically, the Ams distribu¬ 
tion is broad, ranging from Ams = 0 to ~ 30. The median 
ranges from 8 to 13, which nicely brackets the median of the 
ELVIS suite, which is 11. The red, vertical line corresponds 
to Ams = 2 and indicates the number of massive subhaloes 
around the Milky Way, which correspond to the LMC and 
SMC. Based on the 4800 model realizations (which sam¬ 
ple exactly the same host halo masses as the ELVIS suite), 
only 0.8% of the MW-sized haloes have no more than two 
massive subhaloes. The suite-to-suite variance of that per¬ 
centage ranges from 0% to 6%, indicating that a set of 48 
host haloes is insufficient for a meaningful evaluation of the 
TBTF problem. 

3.1.1 Mass and Cosmology Dependence 

The middle panel of Fig. 3 plots the cumulative distribu¬ 
tions of Ams for four 

different host halo masses log[A/o/(h -1 Mq)] = 11.6, 11.8, 
12.0 and 12.2, and for two different cosmologies; ‘WMAP7’ 
and ‘Planck’. The latter has (f2 m ,o, 12 a, o, flb.o, h, as, n s ) = 
(0.3175, 0.6711, 0.0486, 0.6825, 0.8344, 0.9624)’ which are the 
values inferred from the Planck cosmic microwave back¬ 
ground data (Planck Collaboration et al. 2014). Each is com¬ 
puted using a sample of 10,000 model realizations. Note how 
decreasing the mass of the host halo results in a larger frac¬ 
tion of realizations that is consistent with the Milky Way 
in that Ams < 2. Changing the cosmology from ‘WMAP7’ 
to ‘Planck’ causes a slight reduction in the MW-consistent- 
fraction. We stress, though, that these probabilities are very 
sensitive to the exact definition of ‘massive subhalo’. Fol¬ 
lowing Garrison-Kimmel et al. (2014b), in the left and mid¬ 
dle panels of Fig. 3 we have defined massive subhaloes as 
having V aC c > 30 kms -1 and V ma x > 25 kms -1 . In or¬ 
der to gauge the sensitivity to these exact definitions, we 
have repeated the inventory of massive subhaloes chang¬ 
ing the requirement for the present-day maximum circu¬ 
lar velocity to Vmax > 30 kms -1 . This drastically increases 
the MW-consistent-fraction, as is evident from the right- 
hand panel of Fig. 3, which summarizes our results. The 
gray band indicates the probability P(Nms < 2) as func¬ 
tion of halo mass. Solid and dashed lines correspond to the 
WMAP7 and Planck cosmologies, whereas the upper and 
lower bounds correspond to defining massive subhaloes as 
obeying Vmax > 30 kms -1 or 25 kms -1 , respectively. As is 
evident, one can boost the MW-consistent-fraction to well 
over 10% by simply reducing the mass of the MW halo to 
< 10 11 ' 8 A -1 M 0 . 

These results are in qualitative agreement with Wang et 
al. (2012), who, based on an investigation of the Millenium- 
II simulation, argued that the fraction of host haloes having 
three or fewer subhaloes with Emax > 30 kms -1 increases 
from < 5% to ~ 40% as M 0 decreases from 10 1215 ft -1 M 0 
to 10 11 ' 85 h -1 M 0 . Similar results were also reported by 
Vera-Ciro et al. (2013), who, using a semi-analytical model 
of galaxy formation find that if the Milky-Way haloes are 
scaled down to Mo = 10 11 ' 75 h -1 M 0 , the number of satel¬ 
lites brighter than Fornax can be lowered from order of 10 
to 2-5. 
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Figure 4. The distributions of the U max -ratios between the 
and (i — l) th subhalo (in descending order of Umax, where i = 0 
refers to the host halo itself), for i = 1, 2, 3 (from top to bottom). 
The grey distributions are obtained from the 4800 model real¬ 
izations of the ELVIS-size haloes, and the green histograms are 
the ELVIS results. The vertical, red lines mark the corresponding 
MW values, with the red bands indicating the uncertainties due 
to the errors on the V max measurements of the MW, LMC and 
SMC (see Table 1). 

To summarize, according to the massive subhalo formu¬ 
lation the TBTF problem can be significantly alleviated by 
simply lowering the mass of the Milky Way halo to slightly 
below 10 12 /i _1 Mq, which is well within current bounds (e.g., 
Xue et al. 2008; Kafle et al. 2014). However, this alleged 
solution ignores an important observational fact, namely 
that the two massive subhaloes of the Milky Way (i.e., the 
LMC and SMC) actually have fairly large (inferred) val¬ 
ues for Knax- According to van der Marel & Kallivayalil et 
al. (2014) and Kallivayalil et al. (2013), the maximum cir¬ 
cular velocity of the LCM and SMC are 91.7 ± 18.8 kms -1 
and 60.0±5.0kms -1 , respectively. Hence, the MW seems to 
have two subhaloes with Vmax it 60km s -1 . The orange band 
in the right-hand panel of Fig. 3 shows the probability that 
the number, Nmc, of Magellanic-Cloud-like systems is larger 
than or equal to two. The latter are defined as subhaloes 
with Vm ax > 55 kms -1 (upper solid and dashed curves) 
or > 60 kms -1 (lower solid and dashed curves). Clearly, 
the probability P(Nmc > 2) decreases with decreasing halo 
mass. Hence while lowering the mass of the MW halo reduces 
its abundance of massive subhaloes with Unax it 25 kms -1 , 
it also makes it less likely that they have Unax values in 
agreement with the Magellanic clouds. 


3.2 Gap Statistics 

Based on the above findings, the real tension between the 
Milky Way and a simulated MW-size halo is that the third 


most massive satellite has a surprisingly small Umax of ~ 
25km s -1 (irrespective of whether this third system is Draco, 
Ursa Minor, Fornax or the Sgr dwarf) compared to that of its 
second most massive satellite, the SMC, which seems to have 
Vm a x ~ 60kms -1 . Therefore, a more specific formulation of 
the TBTF problem is the existence of a Umax gap between 
the second and the third most massive MW satellites. 

Fig. 4 plots the distributions of the Unax ratio between 
the I th and {i — l) th subhaloes, rank-ordered by Vmax, for 
i = 1,2, and 3. Here the 0 th subhalo corresponds to the host 
halo itself. The model predictions for the 100 mock ELVIS 
suites and the actual ELVIS simulation results are in excel¬ 
lent agreement for all three cases. We take the model dis¬ 
tributions as benchmarks to gauge how (a)typical the MW 
halo is. Using the U ma x values listed in Table 1, the ratio of 
the first subhalo to the host halo, Uma X /Um ax , for the MW 
is 0.54 ± 0.12 (indicated as the vertical, hatched band in 
the upper panel of Fig. 4), well within the theoretical range. 
Similarly, the U max -ratio of the second to first ranked sub¬ 
halo (middle panel) for the MW is 0.65±0.14, again in good 
agreement with the theoretical predictions. Hence, it is clear 
that the subhaloes that host the two Magellanic clouds are 
perfectly common. The TBTF problem is revealed in the 
third panel: assuming Umax = 25 kms -1 for the third MW 
satellite, and the aforementioned Umax measurement for the 
SMC, the V)n a x-ratio for i = 3 for the MW is 0.42 ± 0.04, 
which is small compared to the model prediction whose 95 
percent confidence interval ranges from 0.56 to 0.99. 

Interestingly, the U^ax/V^ ax -distribution predicted by 
the model extends down to values well below that of the 
MW, indicating that it is possible, albeit rare, for dark mat¬ 
ter haloes to reveal a Vmax-gap comparable to that in the 
MW. To make this more quantitative, we construct an en¬ 
semble of 10,000 model realizations for Mo = 10 11,s /i _1 Mg 
in the WMAP7 cosmology and identify realizations with 
Vmfx < 25 kms -1 and Vf^f > 55 kms -1 . Among the 
10,000 realizations, we only find 6 (0.06%) that meet these 
criteria^. The left-hand panel of Fig. 5 compares the cumu¬ 
lative subhalo velocity functions for these six MW-consistent 
cases (solid lines) with that of the actual MW satellites (red 
symbols with error-bars; data taken from Table 1). Note that 
three of the six model realizations mimic the MW extremely 
well down to Umax ~ 15 kms -1 . For smaller Umax all mod¬ 
els predict far more satellites than observed. We emphasize 
though, that the inventory of MW satellites is incomplete at 
the low Vmax end. Tollerud et al. (2008) have demonstrated 
that the luminosity function of MW satellites suffers from in¬ 
completeness for Mv it — 9II. Since there are ten MW satel¬ 
lites brighter than this, we estimate that the MW inventory 
is roughly complete down to Vmax ~ 15kms -1 . Although we 
acknowledge that the rank-order in My is not equal to that 
in Umax, we argue that the discrepancy between model and 
data for Umax < 15 kms -1 is most likely a manifestation of 
incompleteness in the data. 

The right-hand panels of Fig. 5 plot the circular velocity 

If We also repeated this exercise for host haloes with A/q = 
1 q 12 .o /j—i^q, which resulted in only a single candidate. 

I The eight new MW dwarfs reently disovered in the Dark Energy 
Survey (DES Collaboration et al. 2015) are indeed all much fainter 
than this). 
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Figure 5. Rare model realizations with a Enax gap consistent with that observed in the MW. Left-hand panel: grey, solid lines show the 
the cumulative Enax functions for the only six realizations (out of 10,000 haloes with Mo = 10 11 ' 8 h~ 1 Mq, ) that reveal a MW-consistent 
Enax gap. For comparison, the solid circles with error bars show the actual MW data (see Table 1). The median, and the 68, 95, and 
99.7 percentiles of the 10,000 model realizations are indicated in a shaded background with progressively lighter grey. Right-hand panel: 
the rotation curves of the first 12 subhaloes with the highest Enax values in one of those six realizations (the other five cases look very 
similar). The two rotation curves with a thicker line-style indicate the Magellanic-Cloud analogs, while symbols with error bars represent 
the nine brightest MW dSphs (data taken from Wolf et al. 2010). 


curves for the 12 subhaloes with the largest Enax in one of 
our model realizations, randomly chosen from the set of six 
shown in the left-hand panel (results for the other 5 are 
similar). Following Garrison-Kimmel et al. (2014b), these 
circular velocity profiles are computed using the Enax and 
rmax values for each subhalo as predicted by the model (see 
§2.2), and assuming that subhaloes follow an Einasto profile 
(Einasto 1965) 


p(r) = p -2 exp 



( 10 ) 


where r _2 is the radius at which the logarithmic slope of the 
density distribution is equal to —2, p _2 = p(r_ 2 ) and a is 
a shape parameter. In computing our model circular veloc¬ 
ity curves we have adopted a = 0.18, which is the typical 
value for subhaloes in the Aquarius simulation (Springel et 
al. 2008) **. For comparison, the red circles with error-bars 
indicate the circular velocity measurements at the half-light 
radius, E c i r c(?T/ 2 )> of the 9 brightest MW dSphs (Wolf et 
al. 2010). Note how this model predicts two Magellanic cloud 
analogs, with E m ax = 68 kms -1 and 57 kms -1 respectively, 
whose circular velocity curves clearly stand out, with a pro¬ 
nounced gap to the subsequent subhaloes. Note also how 
the circular velocity curves of these subsequent subhaloes 
are statistically consistent with the data for the MW dSphs. 

Hence, based on the gap statistics, we conclude that, 


** Subhaloes in the Aquarius simulation cover the range 0.15 < 
a < 0.30. Note of our results change significantly if we vary a 
over this range. 


within the ACDM paradigm, one can find MW-sized haloes 
with subhalo statistics in excellent agreement with our best 
current understanding of the Milky Way. However, such 
haloes seem to be extremely rare. In order to quantify this 
in more detail, we now investigate the gap statistic in more 
detail, and as function of host halo mass. 


3.2.1 Dependence on Host Halo Mass 

The left-hand panel of Fig. 6 plots the cumulative distribu¬ 
tions of TVgap = N\ — N n . Here Ni and N u are the numbers 
of subhaloes with Enax larger than some lower and upper 
limit, respectively. For the solid lines we set these lower and 
upper limits to be 25 kms -1 and 55 kms -1 , respectively. 
Thus defined, A gap is the number of subhaloes in the range 
25km s -1 < VEax < 55km s -1 . Different colors correspond to 
different host halo masses, as indicated. For each host halo, 
we compute P{< Ngap) using 10,000 model realizations for 
a WMAP7 cosmology. Taking the best-fit values for Enax of 
each MW satellite galaxy, as listed in Table 1, and ignoring 
their uncertainties, the Milky Way has lV gap = 1 (the Sgr 
dSph, for which H m ax = 25.1 kms -1 ). As is evident from 
the left-hand panel of Fig. 6, the probability that a ACDM 
halo has Af gap < 1 is extremely small if Mo >, 10 12 /i -1 Mq, 
but rapidly increases with decreasing Mo- This is quantified 
more clearly by the solid, black line in the right-hand panel 
of Fig. 6 which plots P(N gap < 1) as function of host halo 
mass. We emphasize, though, that these results are very sen¬ 
sitive to how exactly one defines the gap. Setting the lower 
and upper limits of the gap to 30 kms -1 and 60 kms -1 , for 
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example, results in the dashed curves, for which P{N gap < 1) 
is significantly larger. 

The orange, hatched area in the right-hand panel in¬ 
dicates the probability P(N U > 2), where the upper limit 
corresponds to V max = 55 kms' 1 (solid line) or 60 kms -1 
(dashed line), respectively. Note that the probability to host 
two subhaloes consistent with the two Magellanic Clouds 
drops below one percent for Mq < 10 11 ' 4 /i -1 Mg. Hence, to 
be consistent with the MW requires that N gap < 1 and 
N u > 2. The probability that both constraints are satisfied 
is indicated by the red band in the right-hand panel, which 
plots .P(iVgap < 1, N u > 2) as function of halo mass. 

The dashes lines correspond to a different choice for 
the lower and upper limits of 30 kms -1 and 60 kms -1 , re¬ 
spectively. Note how the dashed curves are shifted to the 
left compared to the solid curves, indicating that the typ¬ 
ical iVg a p is significantly smaller for this definition of the 
Knax-gap. Although not shown here, we have verified that 
changing cosmology from WMAP7 to Planck has weak influ¬ 
ence on P(N gap < 1) and P(N U > 2), similar to that shown 
in Fig. 3, and negligible effect on P(N gap < 1, N u > 2). 

The gap statistic reveals a stronger tension than the 
massive-subhalo count. In particular, at best ~ 1% of 
ACDM haloes have a V ma x-gap comparable to that of the 
MW, and that is for a MW host halo mass of Mo — 
10 11,8 /i - 1 M q . This probability drops to below 0.1% for 
Mo < 10 11 ' 2 /i -1 M 0 and Mo > 10 12 ' 2 /i -1 M 0 . These results 
are in excellent agreement with the empirical extrapolations 
of iV-body simulation results by Cautun et al. (2014b). 

3.3 Density of Massive Subhaloes 

The third formulation of TBTF, the one used in the paper 
by Boylan-Kolchin et al. (2011) that introduced the TBTF 
problem, is that the most massive subhaloes in simulations 
are too dense to be consistent with constraints on the MW 
dSphs. Here the internal densities of satellites are usually 
characterized in terms of their r max and 14iax values. 

Measurements of the stellar kinematics of satellite 
galaxies can put accurate constraint on their enclosed mass 
within their half-light radius (Wolf et al. 2010). These in 
turn, constrain a degenerate combination of r ma x and 14,^ 
(e.g., Zentner & Bullock 2003; Boylan-Kolchin et al. 2011). 
This motivated PZ12 to define a density proxy, T, as a linear 
combination of log(r max ) and log(V max ); 

r = 1 + log(0.00141/^ a 2 x /r max ), (11) 

which increases in a direction approximately orthogonal to 
the envelop of the constraint on MW dSphs in the log(r ma x)- 
log(V' max ) plane. Thus defined, T = 1 corresponds to the 20- 
upper bound of the region in the r max -V) nax space that is 
occupied with MW dSphs. A host halo is said to be MW- 
consistent if all its subhaloes obey T < 1, while the presence 
of one or more subhaloes with T > 1 constitutes a manifes¬ 
tation of TBTF. 

The upper left-hand panel of Fig. 7 plots the subhaloes 
in the r max -V max space for the 4800 model realizations of our 
100 mock ELVIS suites, color-coded according to their value 
for Vkcc. The upper right-hand panel shows the same, but 
this time for the actual subhaloes in the ELVIS simulation 
suite of 48 MW-size host haloes. Filled circles indicate the 
median r max as a function of V max , with error-bars indicating 


the 16 and 84 percentiles. Red, dashed lines indicate the 
loci of r = 0, 1 and 2. The model predictions are in good 
agreement with the ELVIS simulation results, in that both 
predict a median r max — Vmax relation corresponding to V > 
1 for Knax 45 kms -1 . Note that all observed MW dwarf 
spheroidals fall in the range 0 < T < 1 (see PZ12). 

Following Garrison-Kimmel et al. (2014b), we fit the 
median R ma x-V max relation with a simple power law 

^max _ , / Vn ax ) P M 

lkpc ll0kms-V ' 1 ’ 

For the model predictions we find (A,p) = (0.77,1.37), in 
good agreement with the ELVIS results, for which ( A,p) = 
(0.73,1.47). A comparison of the best-fit r max -V max relations 
of our model, the ELVIS simulation, and the PZ12 model is 
shown in the lower left-hand panel of Fig. 7. The PZ12 model 
predicts significantly larger scatter and a much shallower 
slope (p = 0.92) compared to both our model and the ELVIS 
simulation. 

The lower right-hand panel of Fig. 7 plots the cumula¬ 
tive distributions of F max = max{Pi}, where the maximum 
is taken over all subhaloes in a single host halo. In the ELVIS 
suite, only one out of the 48 host haloes (corresponding to 
2.1%) has F ma x < 1, and is therefore MW-consistent. This is 
in good agreement with our model predictions, for which we 
find (using a sample of 4800 host haloes) a MW-consistent 
fraction of 2.8%. As in §3.1, we can use the 100 mock ELVIS 
suites to infer the suite-to-suite variance. The median F max 
ranges from 1.1 to 1.3, with a MW-consistent fraction that 
varies from 0% to 15%. This is another demonstration that 
the sample size of the ELVIS suite is insufficient for an ac¬ 
curate assessment of TBTF. 

As is evident from the lower right-hand panel of Fig. 7, 
the PZ12 model predicts a F max -distribution that is much 
broader than what is found in the ELVIS simulation suite or 
predicted by our model. The relatively large difference be¬ 
tween PZ12 and our model is most likely due to two reasons. 
First, PZ12 construct their merger trees using the Somerville 
& Kolatt (1999) algorithm, which is known to overpredict 
merger rates, and result in a halo-to-halo variance of mass 
assembly histories that is too large (e.g., Zhang et al. 2008; 
Jiang & van den Bosch 2014a). Our model, instead, relies 
on the Parkinson, Cole & Helly (2008) algorithm, which 
yields merger trees that are statistically indistinguishable 
from those extracted from numerical simulations. Second, 
PZ12 use a different model for evolving Vnax and r max of 
subhaloes. In particular, whereas we use Eq.(9) to compute 
Vmax and r max , PZ12 simply assume that V max oc m 1 ^ 3 and 
compute r max assuming that dark matter subhaloes follow 
an NFW profile. Based on a number of tests, we conclude 
that this PZ12 model for computing the structural parame¬ 
ters of their subhaloes results in a median r max -V max relation 
that is too shallow, while their merger trees are responsible 
for introducing a scatter in the r max -V max relation that is 
too large, which in turn results in a MW-consistent fraction 
that is too high. 

We caution that the observational constraint on V max 


tt The PZ12 predictions are read off from their published figures, 
which are based on 10,000 realizations of 10 12 /i _ 1 Mq haloes in 
a WMAP7 cosmology. 


© 2014 RAS, MNRAS 000, 1 -19 



12 


Jiang & van den Bosch 



Figure 6. Left-hand panel : the fraction of model realizations with no more than TVcap subhaloes with V m ax E (25,55) kms —1 (solid 
lines) or V mayi E (30, 60)kms —1 (dashed lines), for different host halo masses as indicated (value reflects \og[Mo/(h~ 1 M©)]). Right-hand 
panel, the probabilities of having no more than one subhalo in the V max gap, P(N Q ap < 1), no less than two Magellanic-Cloud analogs, 
P(N u > 2), and both Nq^ p < 1 and N u > 2, PiN q ap < 1, N u > 2), as a function of host halo mass. The solid and dashed curves 
correspond to different threshold V niax values, as indicated in the left-hand panel. 


and r max , expressed as 0 < P < 1, is based on the assump¬ 
tion that dark matter subhaloes have NFW density profiles. 
Vera-Ciro et al. (2013) have shown that if one instead as¬ 
sumes an Einasto profile with a = 0.5, the constraints on 
Vmcix and r max are significantly altered, to the extent that 
even the densest subhaloes in the Aquarius simulations are 
now consistent with the data (i.e., one would no longer infer 
a TBTF problem). However, subhaloes in numerical simu¬ 
lations typically have density profiles that are well fit by an 
Einasto profile, but with a ~ 0.2, for which the constraints 
on Vm a x and r max are very similar to those obtained as¬ 
suming an NFW profile. Hence, the constraint 0 < T < 1 
proposed by PZ12 is still valid, despite the oversimplified 
assumption that subhaloes follow an NFW profile. 

3.3.1 Mass and Cosmology Dependence 

Fig. 8 investigates the dependence of the cumulative F max 
distribution on halo mass and cosmology. The left-hand 
panel plots the results for the WMAP7 cosmology. The 
fraction of realizations with r max < 1 increases from 1% 
at M 0 = 10 12 - 2 /i^Mq to 29% at M 0 = 10 116 /i _1 M 0 . 
Also based on the WMAP7 cosmology, the PZ12 model pre¬ 
dicts the MW-consistent fractions to be 20%, 10% and 10% 
for M 0 = 10 118 10 120 /i^M©, and 10 12 2 r‘M Q 

respectively, significantly higher than our findings of 16%, 
4.5%, and 0.6%. 

The middle panel of Fig. 8 plots the same results but 
now for the Planck cosmology. The MW-consistent frac¬ 
tion increases from 0.1% at Mo = 10 12 ’ 2 /i -1 M© to 9% 
at Mo ~ 10 11,6 /i _1 M@, significantly smaller that for the 
WMAP7 cosmology. Therefore, a relatively small change in 
cosmological parameters seems to have a relatively large im¬ 
pact on the TBTF-statistics in the density formulation. 

Polisensky and Ricotti (2014) argued that the cosmol¬ 
ogy dependence mainly manifests itself as a change in the 
r max of subhaloes. To test this, the right-hand panel of Fig. 8 


plots the r max -V max diagram for the WMAP7 and Planck 
cosmologies, at fixed halo mass of Mo = 10 12 0 /i _1 M 0 . In 
both cases, the r max -Vin ax relations are well fit by Eq. (12) 
with p ~ 1.4. The normalizations of the best-fit relations, 
though, are different, with A = 0.62 for the Planck cos¬ 
mology, and A = 0.74 for the WMAP7 cosmology. This 
indicates that subhaloes in the Planck cosmology are, on 
average, ~ 20 percent denser than in the WMAP7 cosmol¬ 
ogy. The top and side panels of the right-hand panel of Fig. 8 
show the average subhalo F ma x and r max distributions for all 
subhaloes with Knsix > 18kms _1 , respectively. This clearly 
shows that the difference in cosmology predominantly mani¬ 
fests itself as a change in the r max distribution of sublialoes, 
confirming the results of Polisensky and Ricotti (2014). 

The cosmology dependence of subhalo densities arises 
from the cosmology dependence of the host halo assembly 
histories: larger smaller 14 a,o and larger og, as in the 

case of the Planck cosmology compared to the WMAP7 cos¬ 
mology, all result in earlier (average) formation times for 
host haloes of given present-day mass (e.g., van den Bosch 
2002; Giocoli, Tormen & Sheth 2012). Earlier assembly im¬ 
plies that the host halo accreted its sublialoes at earlier 
epochs, when the Universe (and therefore the dark matter 
haloes) was denser. 


4 AN ALTERNATIVE STATISTIC 

In the previous section we have used three different statis¬ 
tics that have been used in the literature to assess the 
TBTF problem. If we adopt a MW host halo mass of 
Mo = 10 12 r'Mg, the inference is that the MW-consistent 
fraction ranges anywhere between ~ 0.1% and ~ 10%, de¬ 
pending on which statistic one uses. Obviously, this raises 
the question which is the more meaningful statistic to use. 
We believe the answer is basically none of the above, and 
the reason is that they either suffer from the ” look-elsewhere 
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Figure 7. Distribution of subhaloes in r max -y max space, for the 4800 model realizations of ELVIS-size haloes (upper, left-hand panel) 
and the 48 ELVIS haloes (upper, right-hand panel). Individual dots are subhaloes, color coded by Vkcc- Circles indicate the median r max 
at given Vmax bins, with error bars indicating the 16 and 84 percentiles. The red dashed lines are the loci of T = 0,1, and 2 [Eq.(11)], with 
0 < r < 1 encompassing the region occupied by MW dSphs. Lower, left-hand panel: comparison of the best-fit r max -V max relations of 
our model (grey), the ELVIS simulation (green), and the PZ12 model (cyan), with solid and dashed lines indicating the median relations 
and the 16 and 84 percentiles respectively. Lower , right-hand panel: corresponding cumulative r max distributions. The red-shaded region 
corresponds to r ma x < 1 and is considered ‘MW consistent’. 

effect” (e.g., Gross & Vitells 2010), and/or disregard certain 
aspects of the data. 

Both the massive subhaloes formulation and the gap 
formulation use statistics that require the identification of 
one or more particular values of Vmax; the massive sub¬ 
haloes formulation considers the number of subhaloes with 
Kn ax above some limit, while the gap statistic is based on 


the number of satellites between two values of V mSLX . These 
values are chosen by the user after carefully examining the 
data on the MW, in an attempt to find a statistic for which 
the MW is least likely. This is a clear example of the look- 
elsewhere effect. For example, if the MW satellite system 
would have revealed a V)n a x-gap between 10 kms -1 and 
30 kms -1 , rather than between 25 kms -1 and 55 knrs -1 , 
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Figure 8. Cumulative T max distributions for different halo masses, in the WMAP7 cosmology (left-hand panel), and the Planck cosmology 
(middle panel). The MW-consistent regime (r m ax < 1) is highlighted in red. The right-hand panel shows the subhalo r max Vmax relations 
for a host halo mass Mo = 10 120 /i — *Mg in the WMAP7 and Planck cosmologies. The circles (WMAP7) and triangles (Planck) correspond 
to a random subsample of model realizations, while the solid and dashed lines indicate the corresponding median relations. The red shaded 
band indicates the region occupied by MW dSphs. The top and side panels plot the P mlx and r max distributions for model subhaloes 
with Vmax > 18 kms -1 . Note that subhaloes are expected to be significantly denser (i.e., smaller r max ) in the Planck cosmology. 


this would have raised a similar concern of being inconsistent 
with ACDM predictions. Yet, such a gap does not manifest 
itself based on the gap statistic used above to assess TBTF. 
Hence, rather than asking what the probability is for a gap 
between 25 kms -1 and 55 kms -1 , one should ponder about 
the probability that a host halo reveals some gap, not nec¬ 
essarily between these two exact values. This is also evident 
from the fact that we have demonstrated that small changes 
in the ‘user-specified’ values that define the gap results in 
large changes in the MW consistent fraction, and thus in the 
inference regarding the severity of TBTF. Ideally, then, one 
should use a statistic that is ‘blind’ in that it does not rely 
on an examination of the data beforehand. 

Another problem with the previous statistics is that 
both the massive subhaloes formulation and the density for¬ 
mulation do not properly account for the Magellanic clouds. 
We believe this to be a serious shortcoming, as the Magel¬ 
lanic clouds, by themselves, put a tight constraint on the 
mass of the Milky Way host halo (e.g., Busha et al. 2011). 

Finally, it is important to realize that no study of TBTF 
to date has properly accounted for the observational errors 
in the y max measurements of the MW satellite galaxies. As 
we demonstrate below, this introduces a huge uncertainty on 
any MW-consistent fraction, and should be properly taken 
into account. 

Based on these considerations, we devise a new statistic 
that is ‘blind’ (i.e., no scale has to be picked upfront), uses 
all data on equal footing, and allows for a straightforward 
treatment of errors in the V max measurements of individ¬ 
ual MW satellites. Consider two rank-ordered distributions, 
Si(xi,X2, ...,xn) and 52( j / i , 3/2, 3/ jv ). In our application, 

Xi and 3 /i are the V m ax values of dark matter subhaloes, 
while 5i and 52 are two different host haloes (i.e., two dif¬ 
ferent model realizations for a host halo of given mass, or 
a model realizations plus the actual F max data for satellite 
galaxies in the MW). Note that 5i and 52 have the same 
number of elements and that Xi +1 > Xi and 3 / 1+1 > 3 /i. We 
now introduce the statistic 


p_ 


(13) 


which is a measure for the difference between the (cumula¬ 
tive) distributions of 5i and 52- In particular, A y\_, \xi — 
3/i | is the absolute value of the area between the cumula¬ 
tive distributions of 5i and 52. We normalize this area by 
jj + 3/i) so that Q is dimensionless, and insensitive 

to an overall shift in x and 3/ (i.e., multiplying Xi and 3/i 
by some factor / leaves Q invariant). Note that Q — 0 if 
5i = S 2 (the distributions are identical), while Q = 1 if 
either 5i or 52 consists solely of null elements (i.e., Xi = 0 
or 3/i = 0 for all i = 1,2, Note that this Q-statistic 

is similar to the Kolmogorov-Smirnov test, which measures 
the maximum value of the absolute difference, dies, between 
two cumulative distributions. However, the KS-test is not 
well suited to characterize differences in the tails of two dis¬ 
tributions; it mainly is sensitive to finding differences in the 
median. We therefore opted to use the statistic Q instead, 
which has equal sensitivity throughout the distributions. In 
adopting Q to assess TBTF, each individual dark matter 
host halo has a corresponding distribution 5, in which the 
elements are the rank-ordered values of V mSLX for the N sub¬ 
haloes with the largest Knax values. 

In order to turn the Q statistic into a probability mea¬ 
sure, we proceed as follows. Given K = 10, 000 model real¬ 
izations, for a given host halo mass, Mo, and a given cosmol¬ 
ogy, we first compute the values Qij for each pair { 5 i, 5 j} 
(with i,j = 1, 2,..., K), where Si is the rank-ordered distri¬ 
bution of the N largest I/ max values for model realization i. 
Next we compute the average 

Qi = j- 'y ' Qij (14) 

j^i 


for each of the 10,000 realizations. Finally, we compute the 
K values of Q mw,i by comparing the I/max distribution of 
the MW to that of each of the 10,000 model realizations, 
which yields 
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Figure 9. Upper panels: The thick, green curves show the cumulative probability distribution, P(> Q ), of the Q statistic defined by 
Eqs. (13)-(14), as obtained from 10,000 model realizations for dark matter haloes in a WMAP7 cosmology. Different panels correspond 
to different host halo masses, as indicated. The solid, red line indicates Qmw? while the red-shaded histogram is the distribution of 
from 10,000 Monte-Carlo realizations of the Vmax-distribution of MW satellites obtained by independently drawing Vmax values from 
their respective error distributions. Blue, downward pointing arrows indicate the Q values for (some of) the model realizations that reveal 
a Vmax-gap similar to the MW [i.e., N(V max > 25 kms -1 ) = N(V m ax > 55 kms -1 ) = 2]. Lower panels: Cumulative Vinax distributions 
for the MW (solid, red line), the median of the 10,000 model realizations (sold, green line), and (some of) the model realizations with 
a Imax-gap similar to the MW. The Q value corresponding to the comparison of the red and green curves is indicated in the top-right 
corner of each panel. 


Qmw 


~K 


K 

y ^ qmw,* 

i= 1 


(15) 


Using the distribution of 10,000 Qi values, we can now com¬ 
pute P{> Qmw), the probability that the ACDM cosmology 
yields host haloes with Q values as large as that of the MW. 
In what follows we refer to this probability as Vmw- 

Finally, we can easily adopt the Q statistic to also ac¬ 
count for observational errors in the U max values of the MW 
satellites. First we construct IVmc = 10,000 Monte-Carlo 
realization of the U ma x-distribution of MW satellites by in¬ 
dependently drawing U max values from their respective er¬ 
ror distributions. For simplicity, we assume that the error- 
distributions for Umax are Gaussian, with mean and stan¬ 
dard deviation equal to the values listed in Table 1. If the 
Unax values quoted in Table 1 has separate upper and lower 
bounds, we assume that the standard deviation is equal to 
the mean of these values (i.e., we ignore any potential skew¬ 


ness in the error distribution). Each Monte-Carlo realization 
results in a set <Sm W , where the prime is used to indicate an 
element of the set of Nmc Monte Carlo realizations. Next, 
for each of these <Sm W we compute Qmw an< A the corre¬ 
sponding Vmw using the same method as described above 
(i.e., by comparing <Sm W to each of the K = 10,000 model 
realizations). The resulting distribution of Pmw indicates 
the uncertainty on Vmw arising from the uncertainties on 
the individual Unax measurements. 

We have performed the above analysis for 8 different 
host halo masses, 

log[Mo/( /i _1 Mq)] = 11.0,11.2, 11.4,..., 12.4, two different 
cosmologies (‘WMAP7’ and ‘Planck’), and 3 different val¬ 
ues for the number of subhaloes in each set, N = 5, 9 and 
20. Fig. 9 shows the results for the WMAP7 cosmology and 
N = 20, for four different values of Mo, as indicated. The 
solid, green curve in the upper panels indicates the cumula¬ 
tive distribution, P(> Q), obtained from the 10,000 model 
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realizations. This P(> Q) is found to be virtually indepen¬ 
dent of host halo mass and cosmology, which is a conse¬ 
quence of the fact that we have normalized Q and that the 
shape of the subhalo V ma x function is largely invariant. The 
shaded histograms indicate the distributions of Qmwi while 
the red, vertical line indicates Qmw- Note how Qmw shifts 
to larger values with increasing host halo mass, and that it is 
always located in the tail of the Q distribution of the model 
predictions, even for halo masses as low as 10 11 ' 4 /? -1 M©. 

The lower panels of Fig. 9 show the cumulative Vnax dis¬ 
tribution, P(> Vmax), of the MW (red histogram), compared 
to the cumulative distribution of the median of the model re¬ 
alizations (green histogram). The latter is obtained by com¬ 
puting the median V ma x of the * th (i = 1 , 2 ,..., 20 ) member 
of all 10,000 model realizations, and roughly corresponds to 
the typical Vnax distribution of the model. The value of the 
Q-statistic corresponding to these two cumulative distribu¬ 
tions is indicated in each panel, and (as expected) is similar 
to Qmw (indicated by the solid, vertical line in the upper 
panels). For a MW host halo mass of Mo ~ 10 12 h -1 M©, 
there is a dramatic discrepancy between model and data: 
Not only does the model predict smaller Vnax for the Magel¬ 
lanic clouds, it also overpredicts Vnax for all other satellites; 
this is a manifestation of the TBTF problem. For a MW 
host halo mass of Mo ~ 10 114 h -1 M©, on the other hand, 
model and data are in remarkably good agreement for the 
satellites with Vmax < 25 kms -1 . However, the model now 
predicts much lower V max values for the Magellanic clouds 
(of the order of 30 — 40 kms -1 ). There is also some tension 
at the low V max end (V ma x ^ 15 kms -1 ), but this reflects 
the onset of the ‘missing satellites’ problem, and may well 
reflect observational incompleteness in the inventory of MW 
satellites. 

The thin, blue histograms show the P(> Vnax) distri¬ 
butions for the model realizations that reveal a Vnax-gap 
similar to the MW (i.e., /V(V m ax > 25 kms -1 ) = N(V max > 
55 kms -1 ) = 2). In the case of Mo = 10 11 ' 8 h - 1 M© 
these are the same 6 models depicted in Fig. 5, while for 
Mo = 10 12 '° /r - 1 M© only 1 (out of 10,000) model realization 
meets these criteria. For Mo = 10 11 ' 4 M© (10 11,6 h - 1 M©) 
we only show a random subset of 10 realizations, from a 
total of 15 (16) out of 10,000 that meet these gap-criteria. 
The Q values corresponding to these model realizations are 
indicated as blue, downward pointing arrows in the upper 
panels. 

In the case of Mo = 10 12 h - 1 M©, it is clear that TBTF 
is not only a problem of an unexpected Vnax-gap between 
the second and third ranked members. After all, the model 
that displays a gap similar to that in the MW (which in 
itself is extremely rare), still is an extremely poor match 
to the Vmax distribution of the other satellite galaxies (at 
least below V ma x ~ 20 kms -1 )! In particular, the Vnax dis¬ 
tribution of MW dwarf spheroidals is much broader than 
what is predicted by the model. Interestingly, in the case 
of Mo = 10 11 ' 4 /i - 1 M©, the models that reveal a MW-like 
Vnax-gap (which has an occurrence rate of ~ 1.5%) basi¬ 
cally are a good match to the entire Vnax distribution of 
dwarf spheroidals, at least down to ~ lOknrs -1 . If it wasn’t 
for the fact that such a low halo mass for the MW is ba¬ 
sically ruled out by the timing argument (see e.g., Phelps , 
Nusser & Desjacques 2013), the space motion of Leo I (e.g., 
Boylan-Kolchin et al. 2013) and the kinematics of blue hor- 



Figure 10. Solid lines show the host halo mass dependence of the 
probability Pmw that the WMAP7 cosmology yields host haloes 
with Q values as large as that of the MW. Open circles and error- 
bars indicate the median and the 18 and 84 percentiles of the 
corresponding distributions, obtained from 10,000 Monte- 

Carlo realizations of the Wax-distribution of MW satellites. Re¬ 
sults are shown for three different values of N, as indicated. 

izontal branch stars (e.g., Xue et al. 2008; Kafle et al. 2014) 
one might be tempted to consider this strong evidence in 
support of a MW halo mass of the order of 10 11 ' 4 /i -1 M©. 

Fig. 10 summarizes the results of our Q-statistic analy¬ 
sis. The solid lines plot the probability Pmw as function of 
host halo mass for three different values of N, as indicated, 
while the open circles and error-bars indicate the median and 
the 18 and 84 percentiles of the corresponding Vm W distri¬ 
butions. These results are all based on the WMAP7 cosmol¬ 
ogy, but the results for the Planck cosmology are virtually 
indistinguishable. If we only focus on the five subhaloes with 
the largest V m ax values, then the inference is that the most 
probable mass for the MW host halo is ~ 10 12 h -1 M©, in 
good agreement with a variety of constraints (e.g., Xue et 
al. 2008; McMillan 2011; Boylan-Kolchin et al. 2013; Kafle 
et al. 2014). Furthermore, the probability that a random 
ACDM dark matter halo of that mass has a V ma x distribu¬ 
tion similar to that of the MW (in terms of the Q-statistic) 
is 6.3jl5 3 3 % (68% CL, when accounting for the errors on the 
individual V ma x measurements of the MW satellites). This 
does not signal concern of a potential problem for the ACDM 
paradigm. 

However, when using N = 20 (i.e., when using the V ma x 
values of all MW satellites in Table 1 down to Willman 1, 
which has Vnax = 8.3l 2 'g kms -1 ), a very different picture 
emerges. In particular, the probability that a host halo of 
10 12 h -1 M© has a Vnax distribution similar to that of the 
MW is now reduced to < 5 x 10 -4 (at 68% CL), while the 
most probable MW mass has dropped to Mo ~ 10 11 ' 4 /i -1 M© 
(where Vmw = 0.6 4 J' 2 % ). As is evident from the lower 
right-hand panel of Fig. 9, the problem is not just a large 
gap between ~ 25 and 55 kms -1 , but rather a problem re¬ 
garding the overall width of the Vnax distribution of the MW 
satellites. Taking into account that our inventory of MW 
satellites may still be incomplete below V ma x ~ 15kms -1 , a 
more robust assessment of TBTF only uses the data of the 
9 highest-Vnax satellite galaxies (LMC, SMC, Sagittarius, 
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Bootes II, Draco, Ursa Minor, Fornax, Sculptor and Leo I). 
In that case, we infer a MW-consistent fraction, for a host 
halo of 10 12 h- 1 Mg, of Vmw = 1.4+11% (68% CL). It re¬ 
mains to be seen how the Vmw for TV = 20 changes as more 
and better data for the population of MW dSphs continues 
to come available. 


5 SUMMARY 

In this paper we have used semi-analytical models for the 
substructure of dark matter haloes to assess the TBTF prob¬ 
lem. We demonstrated that the model accurately reproduces 
the average subhalo mass and velocity functions, as well as 
their halo-to-halo variance, in high resolution TV-body sim¬ 
ulations. We then used the model to construct thousands of 
realizations of MW-size host haloes, which we used to in¬ 
vestigate the TBTF problem with unprecedented statistical 
power. 

Previous studies have formulated TBTF in different 
ways, which we refer to as the massive subhalo formula¬ 
tion, the gap formulation, and the density formulation. We 
have assessed TBTF in all three formulations, and compared 
our results to the ELVIS suite of 48 high-resolution zoom- 
in simulations of MW-sized host haloes (Garrison-Kimmel 
et al. 2014a). Overall, the results from our semi-analytical 
model are in excellent agreement with the ELVIS results. Us¬ 
ing 100 mock ELVIS suites, we demonstrate, though, that 
the suite-to-suite variance is significant, and that a proper 
assessment of TBTF statistics requires simulation suites 
that are an order of magnitude larger than ELVIS. This 
motivates the use of semi-analytical models such as those 
presented here. 

Regarding the three aforementioned formulations, our 
assessment of TBTF is as follows: 

• massive subhalo formulation: according to this for¬ 
mulation, the MW has a deficit of massive subhaloes, defined 
as subhaloes with Umax > 25 kms -1 and Vacc > 30 kms -1 . 
Whereas the MW has only two subhaloes that (most likely) 
meet these criteria, namely the two Magellanic clouds, most 
MW size host haloes have a significantly larger number of 
‘massive subhaloes’, TVms- We find that (TVms) — 7 (8) and 
with P(Nms < 2) ~ 0.1% (0.05%) for a host halo mass 
of Mo = 10 12 Ti -1 M 0 in a WMAP7 (Planck) cosmology. 
However, P(Nms < 2) has a strong mass dependence, in¬ 
creasing to > 10% for Mo < 10 11 ' 8 Ti -1 Mg. In addition, 
it also depends strongly on the exact definition of ‘massive 
subhaloes’. For example, changing the lower limit on U max 
from 25 kms -1 to 30 km s -1 boosts P(Nms < 2) by an order 
of magnitude for Mo = 10 12 Ti -1 Mg. Finally, it is important 
to stress that while lowering the mass of the MW’s halo in¬ 
creases P(Nms < 2), it decreases the probability that such 
a halo hosts two subhaloes comparable to the Magellanic 
clouds. 

• gap formulation: according to this formulation, the 
MW has an unexpectedly large gap in the VT n ax-rank-ordered 
list of its satellite galaxies. In particular, no satellite galaxies 
are known with a Umax between that of the SMC (U max = 
60 ± 5 kms -1 ) and that of the Sagittarius dSph (U ma x = 
25.1±1.5km s -1 ). If we define ‘MW-consistent’ as having two 
subhaloes with U ma x > 55 kms -1 and at most one subhalo 
with 25 < Umax < 55 kms -1 , then we find that the MW 


consistent fraction peaks at ~ 0.6% around a host halo mass 
of Mo — 10 11 ' 8 /i -1 M 0 . For Mo > 10 12 Ti -1 Mg this fraction 
is found to be less than 0.1%, with little dependence on 
cosmology. As for the massive subhalo count, though, these 
MW consistent fractions are extremely sensitive to the exact 
definition of the gap; for example, changing the U ma x-values 
of the gap by a mere 5km s -1 can change the MW-consistent 
fraction by an order of magnitude. 

• density formulation: according to this formulation, 
the densities of the (more massive) subhaloes are too high 
compared to those of the subhaloes hosting MW dSphs. In 
particular, PZ12 introduced a density parameter, T, defined 
by Eq. (11), and argued that whereas all MW dSphs have 
0 < r < 1, the majority of MW-size host haloes have at 
least one sublialo with T > 1 (such a subhalo is consid¬ 
ered too big, or rather dense, to fail). We find that, in a 
WMAP7 cosmology, only ~ 4.5% (15.8%) of host haloes 
with Mo = 10 12 Ti -1 Mg (10 11 ' 8 Ti -1 Mg) are MW-consistent 
in that r max < 1. Note that PZ12, using a semi-analytical 
model similar to ours, claimed that these MW-consistent 
fractions are significantly higher, at 10% and 40%, respec¬ 
tively. As discussed in the text, we believe that their re¬ 
sults are hampered by the use of inaccurate halo merger 
trees, and the oversimplified assumption that subhaloes have 
NFW density profiles. Finally, we emphasize that of all three 
TBTF formulations, the density formulation is the one most 
sensitive to cosmology. In particular, we find that, in a 
Planck cosmology, the MW-consistent fraction of host haloes 
with Mo = 10 12 Ti -1 Mq is only 0.6%, rather than 4.5% in a 
WMAP7 cosmology. 

An important, and troubling, downside with all three 
formulations above is that they suffer from the look- 
elsewhere effect and/or disregard certain aspects of the data 
on the MW satellite population. In an attempt to rem¬ 
edy these shortcomings, we have devised a new statistic 
which is ‘blind’, in that it doesn’t require the pre-selection 
of some particular scale (i.e., a particular Unax, Vacc and/or 
F value), treats all data on equal footing, and also allows 
for a straightforward treatment of errors in the U ma x mea¬ 
surements of individual MW satellites. Similar in spirit to 
the Kolmogorov- Smirnov (KS) test, this Q-statistic, defined 
by Eq. (13), compares two cumulative distributions. It pro¬ 
vides a (conveniently normalized) measure for the absolute 
value of the area between two cumulative distributions, and 
has the advantage over the KS-test that it has equal sensi¬ 
tivity throughout the distributions. As a consequence, it is 
equally sensitive to the offset between two distributions as 
to a difference in the spread. 

Using this Q-statistic to compare the U max distribu¬ 
tion of the 9 MW satellites with the largest Unax measure¬ 
ments to that of dark matter subhaloes in our model real¬ 
izations, we infer a MW-consistent fraction, for a host halo 
of 10 12 fc -1 Mg, of Vmw = 1-4 7 8 ' 8 %, where the upper and 
lower bound reflect uncertainties due to the errors on the 
individual U ma x errors. This is similar to the MW-consistent 
fraction inferred in the density formulation, but significantly 
larger than what is inferred from the gap-statistic. However, 
the latter is severely hampered by the look-elsewhere effect, 
and does not reflect a proper assessment of the TBTF sever¬ 
ity. 

To conclude, it is difficult to express the severity of 
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TBTF in a single number. In general, TBTF is (slightly) 
more problematic in a Planck cosmology, compared to a 
WMAP7 cosmology, and is minimized if the virial mass of 
the MW’s host halo falls in the range 3 - 6 x 10 11 /i -1 Mq. 
Such a low MW mass, however, is basically ruled out by 
a variety of independent constraints (e.g., Xue et al. 2008; 
McMillan et al. 2011; Boylan-Kolchin et al. 2013; Phelps et 
al. 2013; Kafle et al. 2014). Assuming instead a host halo 
mass of 10 12 and using only data on the 9 known 

satellite galaxies with Knax > 15 kms -1 , both the density 
formulation and the Q-statistic suggest a MW-consistent 
fraction of the order of a few percent, something which 
we do not consider particularly challenging for the ACDM 
paradigm. However, if it turns out that the inventory of 
MW satellite galaxies is complete to ~ 8 kms -1 (i.e., fu¬ 
ture surveys uncover no new MW satellite galaxies with 
Vmax >, 8 kms -1 ), then it is clear that the spread in 14,^ 
values for the 20 highest Vmax-ranked satellites is utterly in¬ 
consistent with ACDM predictions. In that case, we either 
(i) have systematically underestimated the V)nax values (by 
about a factor of two) of all MW dwarf spheroidals, (ii) the 
process of galaxy formation significantly lowers the central 
densities (and hence V ma ,J) of the subhaloes hosting dwarf 
spheroidals, or (iii) the ACDM paradigm is actually falsified. 
In this respect, it remains to be seen whether baryonic ef¬ 
fects, such as those discussed in Zolotov et al. (2012), Brooks 
& Zolotov (2014), and Arraki et al. (2014), can give rise to 
satellite populations in hydro-dynamical simulations with 
Vmax distributions in better agreement with observations. 
We hope that the Q-statistic introduced here will prove use¬ 
ful to compare such simulations to the data. 
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